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Abstract 

We determine numerically the full complex Lyapunov spectrum of SU(2) Yang-Mills fields 
on a 3-dimensional lattice from the classical chaotic dynamics using eigenvalues of the mono- 
dromy matrix. The microcanonical equation of state is determined as the entropy - energy 
,_h \ relation utilizing the Kolmogorov-Sinai entropy extrapolated to the large size limit. 

Hi ' 

C^l ■ Introduction 

Knowledge of the microscopic mechanism responsible for the local equilibration of energy 

qq ' and momentum in the nuclear collision is fundamental to understand the possible emergence 

of a quark gluon plasma in relativistic heavy-ion collisions. It is evenly essential to study 

the equation of state of this system (in and out of equilibrium) in order to be able to make 

predictions or interpret experimental data from the perspective of quark - gluon matter. 

Many of the approaches to the equation of state assume thermal equilibrium, since the real 

■^> ' - time evolution of such a plasma in the full quantum field theory cannot yet be followed in a 

non-perturbative manner by presently practized numerical methods. Perturbative QCD helps 

to get hints from the high-energy parton dynamics, but the non-perturbative lattice gauge 

C^h, theory is difficult to extend to non-equilibrium phenomena. 

^ '■ 

Fortunately, the situation is eased by findings which point to the possibility that the quark 

- gluon plasma can be dynamically studied in a semiclassical approximation. Hard thermal 
k> , loop resummation has been shown to cope with the semiclassical transport (linear response) 

theory approach [EL 0] . In fact earlier numerical studies of the classical Hamiltonian dynamics of 
Yang-Mills systems on three - dimensional lattices revealed an intriguing coincidence between 
the maximal Lyapunov exponent and twice the gluon damping rate B. In the background 
of that argumentation equipartition of the energy has been conjectured as a result of the 
classical chaotic dynamics on the lattice. This conjecture has been supported by the pseudo 

- Boltzmannian distribution of one-plaquette energies (the local magnetic energy density) in 
long-evolved classical configurations. 

The correspondence between this average energy and the Kolmogorov-Sinai entropy has been 
first investigated in 0] for pure SU(2) Yang-Mills systems. Those times, however, the Lyapunov 
spectrum could be obtained only for relatively small systems (N=2,3) with the rescaling method. 
This way also only the positive real exponents could have been calculated. The question of the 
thermodynamical limit (N — > oo) and the investigation of a possible dependence on the initial 
configuration (which is usually done by using improved statistics) remained open. (Citation: 
"Future work is still required here" [||, page 6].) 

In the present article we study the ergodization of SU(2) lattice gauge theory due to its 
classical chaotic dynamics by using different methods as before. We also consider larger lattices 
(N = 2, 3, 4, 5, 6) for obtaining the entropy and extrapolate to the large N limit. In our approach 
the full complex spectrum of Lyapunov exponents will be obtained using the monodromy matrix 



(linear stability matrix along the path determined by the classical equation of motion) in the 
space of evolving field configurations on the lattice. The complex spectrum offers an insight 
also into the periodic (oscillatory) behavior of field fluctuations; a basic experience for learning 
about the corresponding quantum theory. 

The question of ergodization is addressed via the Kolmogorov-Sinai entropy, which is cal- 
culated up to linear lattice sizes of N — 6 (using 24N 3 = 5184 dimensional phase-space) and 
extrapolated to N — ► oo. We extrapolate both for long times, - in order to get close to the 
Lyapunov exponents, - and for large systems - in order to approach the thermodynamical limit. 
Ensembles of randomly chosen (i.e. chosen according to the Haar measure in the group vari- 
ables) initial configurations with a given (slightly fluctuating) total energy are averaged over, 
too. 

This article is constructed as follows. First we review basic definitions dealing with the 
chaotic behavior of extended systems and then the Hamiltonian treatment of classical lattice 
gauge theory. Finally we present results on full complex Lyapunov spectra, scaling properties, 
extrapolations and the entropy - energy relation (microcanonical equation of state) for SU(2) 
lattice gauge theory. 



Measuring Chaos 

Instead of the classical determination of the Lyapunov exponent, which has been calculated 
by the rescaling method in ||, |[ g, 0, g] , in this article we use the monodromy matrix approach. 
The Lyapunov spectrum Li is expressed in terms of the monodromy matrix's eigenvalues A$: 

U = hm £^ * = 1..../ (1) 

where the A^(i)-s are the solutions of the characteristic equation 

det(Ai(t)l-M(t)) = (2) 

at a given time t. Here M is the linear stability matrix, / is the number of degrees of freedom. 

We consider conservative dynamics fulfilling Liouville's theorem. 

/ 

E^ = ° ( 3 ) 

i=0 

In the numerical calculation we use the discrete definition of the Lyapunov spectrum: 

n 

L l = (Ai} W = -EAife-i) i = W, (4) 

U 3=1 

where tj are subsequent times along an evolutionary path of the gauge field configurations. Nev- 
ertheless a few initial configurations has been omitted from the average, although the starting 
configurations has been produced randomly. We need to extrapolate the quantities L i to the 
long time limit n — > oo, with fixed timestep. It is expected to converge to Lyapunov exponent 
Li occuring in eq. O) for a non-compact configuration space. As we shall see by discussing 
the results of numerical simulations gauge group systems live in a compact configuration space, 
therefore this limit is not entirely safe in our case. The rescaling method observes pairs of 
gauge field configurations which are close in the phase space, so it has a control on the distance 
of these two and can monitor this quantity not growing close to or over the limiting distance 



given by the compact size of the phase space. It is achieved by doing frequent rescalings. The 
monodromy matrix method on the other hand follows only one gauge field evolution, in this 
case the short and long time behavior can be different. 

The Kolmogorov-Sinai entropy is obtained by using Perin's formula: 
h KS = Y,LS{L i ) (5) 

i 

using the theta function Q(x) being one for positive arguments and zero otherwide. The di- 
mension of h KS yet is a rate (1/time), so for estimating the entropy we shall use the normalized 
quantity: 

h KS 
S= Re(L )Nf (6) 

The redefinition is done in a way typical for extensive quantities, here the entropy density, on 
an N x N x N lattice. 



Lattice Yang-Mills field systems 



Our discussion is based on the Hamiltonian formulation of the classical lattice SU(2) gauge 
theory. It is governed by the Hamiltonian M : 



h = J2 (\ (u*,i, u*,i) + ( 1 -\ ( u ** V ^n ) 



(7) 



where U Xj i is the SU(2) group element lying on an oriented link of the lattice starting at the site 
x, pointing to the i direction. In the numerical realization it is represented by real quaternions. 
The notation < A. B > belongs to the normalized trace of the product of two SU(2) group 
elements : 

<A,B >= -ir(AB _1 ). (8) 

The V x .i complement matrix is constructed from triple products of link variables, which com- 
plete the considered link (x,i) to an elementary plaquette. 

V*,i = \Y,U i Uf 1 U£ 1 (9) 

Finally the overdots denote the derivative with respect to the scaled time, t/a, as well as the 
Hamiltonian H stands for the scaled energy g 2 aE. 



The Hamiltonian equations of motions are given by (suppressing the link indices in the 
rotation): 

U = P, 

P = V-(U,V)U-{P,P)U. (10) 

These equations of motion conserve unitarity < U, U >= 1 , orthogonality < U,P >= 0, and 
Gauss' law. The monodromy matrix is defined as: 

dii du 

dU dP 

M=\ | (11) 

dP dP 
dU dP 



Here the different partial derivatives are obtained by comparing the above equations of motion 
for the configurations (U, P) and (U + SU, P + SP). We arrive at 



dU a 
811* 
8U a 



= 0, 
= S ab , 



Qpb 

BP a 8V a f <W C \ 

^TK = -rrff" U c —)u«-V b U«-(U c VC + P c n5»\ 



dU b dU b \ c dU b J 

f)pa 

— b = -2P b U". (12) 

These expressions provide information about the stability of trajectories in the neighbourhood 
of any point of an orbit in the (U,P) phase space. A small perturbation, (SU, SP) evolves in 
time governed by the monodromy matrix M. The eigenvalues of this matrix can be classified 
as follows: for real and positive eigeinvalues neighbouring trajectories apart exponentially and 
the motion is unstable. In the limit of large time we obtain the Lyapunov exponents from these 
eigenvalues. The imaginary parts of the complex eigenvalues describe oscillatory frequencies of 
perturbations. 



In contrast to the rescaling method, where only the positive real part of the Lyapunov 
exponent can be measured, the study of the spectrum of the monodromy matrix makes these 
complex eigenvalues numerically available for the first time. The total number of degrees of 
freedom in the numerical simulation is given by / = 4*3* N 3 — 12 *N 3 , using four- real element 
quaternions for SU(2) group elements (so the phase space is 2/ = 24A 3 dimensional). Due to 
the constrains < U,U >= 1, however, the physically relevant number of degrees of freedom 
is the same as in earlier simulations, which used E a , and tr(C/r a ) as independent degrees of 
freedom. 



The spectrum of the stability matrix 



The eigenvalue spectrum Aj of the monodromy matrix consists of a real subset and a complex 
subset on the complex plane (cf. Figs. 1 and 2). Although this 2/ x 2/ sized matrix is sparse, 
since each group element exerts a force onto close links (those sharing a common plaquette) 
only, the usual sparse matrix methods are not applicable, because they determine just a few 
leading eigenvalues. In order to determine all complex eigenvalues with an acceptable precision, 
only 'brute force' methods with a memory need of 0(f 2 ) and computational time of C(/ 3 ) can 
be used. The maximal system size, our computational resources allow us to use, belongs to 
N = 6, meaning a 2/ = 5904 dimensional phase space. 



Figs.l and 2 display the complex eigenvalues for several configurations during the evolution 
along a few trajectories. At high energy per degree of freedom, g 2 aE — 0.8, (near to the 
saturation value 1) the region covered by the eigenvalues in the complex plane shows a mirrorred 
"butterfly" shape, while at low energy, g 2 aE = 0.1, rather two "ovals". In all cases the covered 
region is symmetric both to the real and imaginary axes: the former property is due to the fact 
that the equations of motion are real, the latter is due to the Hamiltonian is conservative (time 
independent). It is interesting to observe that about one sixth of the eigenvalues are purely 
real, not allowing for oscillations or wave-propagation in the fluctuations. In neither case shows 
the eigenvalue region any resemblence of Wigner semicircles, like those found in studies of the 
eigenvalues of the Dirac operator in random gauge field background |l0| , and usually regarded 



as a signature of quantum chaos. Here we deal with a bosonic system at comparatively high 
excitations and the chaotic behavior is quite classical. 



According to the imaginary part of the eigenvalues the pattern separates to three islands 
at all energies. An important qualitative feature of these patterns is the gap between zero and 
non-zero imaginary parts: it behaves as a dynamically developed infrared cut-off ( "gluon mass" ) 
for oscillations in small perturbations of the classical equations of motion. Also a number of 
zero-frequency modes occur, this is connected to symmetry transformations commuting (in the 
Poisson bracket sense) with the Hamiltonian, such as time independent gauge transformations. 
The gap seems to be reduced and eventually disappear at high energy. Do these differences 
between low and high energy establish a two-phase picture of the classical lattice SU(2) system? 
In order to answer this question a study of the equation of state is necessary. 



Scaling and Equation of State 



We aim to obtain the equation of state from the dynamical simulation. For this purpose we 
first study the finite size scaling and extrapolate to infinite (1/N — 0) lattices. Then we study 
the energy dependence of the maximal Lyapunov exponent as well as the Kolmogorov- Sinai 
entropy. The latter leads eventually to the equation of state as the entropy - energy relation, 
S(E), in the thermodynamical limit of infinite volume. 



Fig. 3 shows the real part of the Lyapunov spectrum extrapolated to 1/N — > from data 
taken at N = 2,3,4,5 and 6 at high energy (g 2 aE = 0.8). The overall pattern is similar 
to that obtained earlier for smaller systems (N = 2,3), just the number of purely imaginary 
eigenvalues (for those ReLi — 0) is greater due to our use of more variables (and, of course, 
more constraints) . The structure of the ordered real part of the Lyapunov spectrum is similar 
at all energies, but the maximal point, clLq, scales with the energy, g 2 aE. 



Fig. 4 displays the fluctuation of the maximal real part Lyapunov exponent, Lq, originating 
in different, randomly chosen initial field configurations. This levels at a few per cent. The 
energy scaling of the maximal Lyapunov exponent has been debated in the past and has been 
found linear in the long-time limit. Doubts rose for low energies [O, H2] , stating that the correct 
scaling here would be L$ ~ E 1 ' A . More elaborated studies with the rescaling - aparting method 
showed then a tendency back to the linear, L$ ~ E scaling also at low energies for long enough 
time observations O, [Ul Hi. Our data agree with the linear scaling with a coefficient of 0.5 
for a middle to short time evolution, as it can be seen in Fig. 5. Surprisingly, and not yet fully 
understood, the very long time behavior of the maximal real part of the monodromy matrix 
eigenvalues show a definitely sublinear scaling with the energy per degree of freedom. The best 
fit is, however, not like E 1 / 4 , but rather a logarithm (cf. Fig. 6). We suspect that following a 
trajectory too long makes the observed eigenvalues feel the compactness of the configuration 
space - which is an artefact of the lattice field theory Hamiltonian. In the following discussion 
we refer to data with linear energy scaling only. 



Fig. 7 displays the extrapolation of the maximal Lyapunov exponents to the thermodynam- 
ical limit at different energies. The correspondence proved to be almost linear by assuming 
an 

Lo ~ _L „ iV-3/2 

v 7 / 



scaling with the finite size. This corresponds to sampling ergodic states 15 . In particular at 
high energies the extrapolated Lq is higher than any actually obtained value in simulations at 
the finite N. 



Finally we obtained the Kolmogorov-Sinai entropy from infinite-size extrapolated, initial 
state and evolution averaged eigenvalue spectra as a function of the scaled energy. For a 
(nearly) ideal gas on expects an S ~ log E relation, indeed this is a good approximation of our 
data (cf. Fig. 8). We obtain 

(S) = (0.5084 ± 0.023) log( 5 2 a.E) + (2.3334 ± 0.0452). (13) 

The indicated errors after the ± signs are of statistical nature. This best fit belongs to an 
inverse temperature 



1 

T ~~ 


d(S) , 
dE 


^ °- 5 

~ ~E 




describing equipartion 


with 


E = 


■¥■ 







(14) 



(15) 

Based upon our calculations even a two-phase structure (or a crossover) cannot be excluded with 
absolute safety: at mid energies a depletion is hinted in the data. The S{E) curve would show a 
first order two-phase structure by having a break somewhere. For a possible relevance to lattice 
SU(2) systems see the paper [jig . In our case, starting from the chaotic dynamics of classical 
lattice Hamiltonians, a refined study is needed in the transition region g 2 aE = 0.2, . . .0.5 in 
order to settle this question. 
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FIGURES 



SU(2) complex Lyapunov spectrum 




Figure 1: Full complex eigenvalue spectrum of monodromy matrix: high energy. 



SU(2) complex Lyapunov spectrum 




Figure 2: Full eigenvalue spectrum of the monodromy matrix: low energy. 
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Figure 3: The real part of the Lyapunov spectra in the case of N — > oo 







N=6, E=0.964, 10 different seeds 
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Figure 4: Maximal Lyapunov exponent obtained starting from different random configurations 
on an TV = 6 lattice at scaled energy g 2 aE = 0.964. 
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Lyapunov scaling, N=6 
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Figure 5: Scaling of short time Lyapunov exponents with energy (before reaching the satu- 
ration in the distance of two initially adjacent configurations due to the compactness of the 
configuration space). 
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Lyapunov scaling, N=6 



fit: f(x) = 0.199 * log(1 +X/0.051) 
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Figure 6: Scaling of Lyapunov exponents with energy for long time evolution of the monodromy 
matrix (after reaching distance saturation in the compact configuration space). 



0.8 
0.7 
0.6 

0.5 

X 

J 0.4 
0.3 
0.2 

0.1 








Finite size scaling 
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Figure 7: Finite size scaling of the maximal Lyapunov exponent. The straight line fits belong 
to different scaled energies (marked in the plot). 
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SU(2) Yang-Mills EOS 
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Figure 8: The normalized Kolmogorov-Sinai entropy. The best fit belongs to \kT energy per 
degree of freedom, but a two phase structure also cannot be excluded (depletion hint). 
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